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ABSTRACT 

We used Herschel Hi-GAL survey data to determine whether massive young stellar 
objects (MYSOs) are resolved at 70/im and to study their envelope density distri¬ 
bution. Our analysis of three relatively isolated sources in the l = 30° and l = 59° 
Galactic fields show that the objects are partially resolved at 70/im. The Herschel 
Hi-GAL survey data have a high scan velocity which makes unresolved and partially 
resolved sources appear elongated in the 70/im images. We analysed the two scan direc¬ 
tions separately and examine the intensity profile perpendicular to the scan direction. 
Spherically symmetric radiative transfer models with a power law density distribution 
were used to study the circumstellar matter distribution. Single dish sub-mm data 
were also included to study how different spatial information affects the fitted density 
distribution. The density distribution which best fits both the 70/im intensity profile 
and SED has an average index of ~ 0.5. This index is shallower than expected and is 
probably due to the dust emission from bipolar outflow cavity walls not accounted for 
in the spherical models. We conclude that 2D axisymmetric models and Herschel im¬ 
ages at low scan speeds are needed to better constrain the matter distribution around 
MYSOs. 

Key words: circumstellar matter - infrared: stars. 


1 INTRODUCTION 

The formation of massive stars presents many challenges due 
to the competing and interlinked roles of gravity, magnetic 
fields and radiation. It is becoming clear through numeri¬ 
cal simulations that material can continue to accrete on to 
a luminous, massive forming star via an accretion disc de¬ 
spite the strong radiation pressure on dust (Krumholz et al. 
2009; Kuiper et al. 2010). Bipolar outflows appear to be a 
ubiquitous ingredient in the star formation process driven 
by magnetic forces (Banerjee & Pudritz 2007) which also 
helps relieve the extreme radiation pressure (Cunningham 
et al. 2011). 

These competing infall and outflow processes shape the 
circumstellar matter distribution around massive forming 
stars. The corollary of this is that a detailed mapping of the 
circumstellar matter distribution can be used to constrain 
models of massive star formation. One way to probe the cir- 
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cumstellar matter close to the protostar is to use the emis¬ 
sion from the heated dust. This has some advantages over 
using molecular line emission for which complex chemical 
and excitation effects have to be taken into account before 
the total gas density distribution can be recovered. There 
are also disadvantages with using warm dust emission as it 
does not convey any kinematic information and becomes op¬ 
tically thick for A < 100/im. However, a full understanding 
of the dust emission also yields the temperature distribution 
of the material which is an important input back into the 
molecular line diagnostic process. 

Different IR wavelengths will probe regions at differ¬ 
ent distances from the central accreting source due to the 
temperature gradients. Typical temperature gradients vary 
from T oc r -3 / 4 in the optically thick part to T oc r -1 / 2 
in optically thin regions (eg. Ivezic & Elitzur 1997). Taking 
the over-simplified approach of the Wien Displacement law 
to locate where most of the emission arises from, we then 
see that the size of the emitting region r oc A 4 ^ 3 or r oc A 2 in 
optically thick and thin regions respectively. In the immedi¬ 
ate environment of a protostar we will have optically thick 
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conditions in the dense mid-plane regions whilst the bipolar 
outflow cavities will be mostly optically thin. 

High resolution studies of the warm dust emission 
around massive forming stars have included the 8-13/rm 
interferometric studies by de Wit et al. (2007;2010;2011). 
These 40 milliarcsecond resolution studies probe size scales 
of about 100 AU for typical distances of nearby massive 
young stellar objects (MYSOs) which approaches the size 
of the dust sublimation radius of about 25 AU. The mid-IR 
visibilities are mostly matched by 2D axisymmetric radia¬ 
tive transfer models where most of the emission arises from 
the warm dust along cavity walls. A compact element inside 
the dust sublimation radius such as an accretion disc may 
be needed to explain the rising 8/xm visibilities. 

Since the dusty bipolar cavity walls are directly illumi¬ 
nated by the central star then this is mostly optically thin 
emission. We would therefore expect the size of the emitting 
region to scale as A 2 . For single-dish observations then this 
size scale is getting larger faster than the diffraction-limited 
resolution is degrading. Hence, for a given diameter tele¬ 
scope it is better to use the longest wavelength possible. For 
the thermal-IR regime from the ground this is 24.5/xm at the 
far end of the Q band. Scaling from the 60 milliarcsecond 
size for W33A at 13/xm (de Wit et al. 2007) this would yield 
a size of 0.2" at 24.5/xm which would be partially resolved by 
the 0.6" diffraction limit for 8 m telescopes. This is what was 
found by de Wit et al. (2009) and Wheelwright et al. (2012) 
who partially resolved a sample of massive YSOs. de Wit 
et al. (2009) modelled the extended emission with spherical 
models and concluded that the density distribution needed 
to be n oc r -1 . This was interpreted as requiring some ro¬ 
tational support on these scales as it is shallower than the 
n oc r -1 ' 5 expected for free-falling material. Wheelwright 
et al. (2012) again used 2D axisymmetric radiative transfer 
models to show that the 20/im emission is also dominated 
by the warm envelope dust along the cavity walls. This is 
seen explicitly in some edge-on systems (e.g. De Buizer 2005; 
2006). 

The recent release of Herschel space telescope data pro¬ 
vides the highest resolution images to date at far-IR wave¬ 
lengths. If we extend the scaling above to the shortest wave¬ 
length of Herschel of 70/xm then we would expect a size of 
the emitting region of about 2". This is comparable to the 
5" diffraction-limited beam of Herschel at 70/xm. Hence, we 
expect to recover further spatial information on the some¬ 
what cooler dust located further from the central object, but 
not so far as to be more influenced by the general molecular 
core environment and ambient radiation held as is likely at 
the longer Herschel wavebands. 

The most extensive observations of massive forming 
stars with Herschel comes from the Hi-GAL survey of the 
Galactic Plane (Molinari et al. 2010a) where most massive 
young stars are located. Here we examine these data from 
the first two Science Demonstration Phase fields observed as 
part of Hi-GAL (Molinari et al. 2010a) to determine whether 
massive YSOs are indeed extended at 70/im. 

We define massive YSOs as deeply embedded, lumi¬ 
nous (L>3000L©), mid-IR-bright, point sources that are 
not ionizing their surroundings to form an ultra-compact H 
II region (UCHII). The lifetime of this phase is about 10 5 
years (Davies et al. 2011; Mottram et al. 20111a). Davies 
et al. (2011) show that the luminosity function of MYSOs 


and UCHII regions is consistent with the MYSOs becoming 
swollen due to high accretion rates as predicted by the mod¬ 
els of Hosokawa V Omukai (2009). This means their effective 
temperatures are too low to ionise their surroundings until 
either they stop accreting at high rates or grow to greater 
than about 3OM0, when they contract rapidly down to their 
main sequence radius. 

The sample of MYSOs we use comes from the Red MSX 
Survey (RMS) (Lumsden et al. 2013). Starting from an ini¬ 
tial colour selection of mid-IR bright sources (Lumsden et al. 
2002) from the MSX satellite point source catalogue (Price 
et al. 2001) we have followed these up with radio contin¬ 
uum observations to distinguish UCHII regions and dusty 
planetary nebulae (Urquhart et al. 2009); determined kine¬ 
matic distances from 13 CO observations and H I absorption 
(Urquhart et al. 2008; 2012); and determined luminosities 
from far-IR fluxes (Mottram et al. 2010; 2011b). All MYSOs 
discussed in this paper are undetected at 5 GHz at the 1 mJy 
level. 

In this paper we examine the Herschel imaging of the 
RMS MYSOs in two lots of 2 square degree regions of the 
Galactic plane. The peculiarities of the Hi-GAL 70/xm point 
spread function (PSF) are discussed in section [2] and the 
70/xm imaging of the RMS MYSOs in these regions is de¬ 
scribed in section [3j Radiative transfer modelling of the 
70/mi emission is presented in section [H and the results dis¬ 
cussed in section [5] Conclusions are drawn in section [6] 


2 HI-GAL 70/im POINT SPREAD FUNCTION 

The Hi-GAL survey was carried out in a specially developed 
parallel mode whereby the PACS and SPIRE instruments 
simultaneously scan the sky at a fast rate of 60" per second 
(Molinari et al. 2010). This causes the nominally circular 
PACS 70/xm beam with FWHM of 5.3" to be smeared out in 
the scan direction with a resultant size of 5.8" x 12.1" (Lutz 
2012). A high signal-to-noise representation of the PACS 
parallel mode 70/xm PSF is shown in figure [I] which shows 
an image of the asteroid Vesta observed in similar conditions 
as Hi-GAL data (see Lutz 2012). The first Airy ring can also 
be seen smeared out along the scan direction (PA=42.5°) 
and a dark spot is seen at one end of the scan direction. 

The two Science Demonstration Phase (hereafter SDP) 
Hi-GAL fields were visually inspected to search for objects 
that might be suitable PSF objects at 70/im, i.e. bright, 
unresolved and isolated. Such objects are rare and we only 
found two in each of the SDP fields and their details are 
listed in table [T| These appear to be all AGB or post-AGB 
stars and they were all located away from the mid-plane 
consistent with an evolved, intermediate mass population. 
Such stars are losing mass with dusty winds which makes 
them suitable for IR PSF stars. However, as de Wit et al. 
(2009) found these stars can also be extended and so have 
to be treated with caution when using as PSF objects. 

An image of the PSF object V1362 Aql from the 1=30° 
SDP field is shown in figure [2] This is from the so-called 
naive map constructed from the Hi-GAL data which adds 
together the two different scan directions referred to as nom¬ 
inal and orthogonal. Unfortunately this results in a complex 
PSF which is basically a cross-shape. The naive maps we 
used had not had any astrometric corrections applied, which 
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Figure 1 . Image of Vesta taken in parallel mode with a scan 
speed of 60 // per second and array-to-map inclination angle of 
a = 42.5° showing the details of the point spread function. In 
order to use a logarithmic stretch, a constant value of 10 was 
added to the normalised fluxes. 


Table 1. Parameters of the PSF objects found within the two 
Hi-GAL SDP fields at 70/un. 


Name 

Nature 

RA 

Dec 

/70^m 



(J2000) 

(Jy) 

V1362 Aql 

Mira 

18:48:41.9 

-02:50:28 

66 

IRAS 18491-0207 

PAGB 

18:51:46.2 

-02:04:12 

80 

IRAS 19374+2359 

PAGB 

19:39:35.5 

+24:06:27 

29 

IRAS 19348+2229 

? 

19:36:59.8 

+22:36:08 

32 


also resulted in an offset cross-shape. Such a complex PSF 
makes the search for extended emission very difficult and 
certainly precludes the use of azimuthally averaged radial 
profiles as we used for ground-based 24/xm imaging (de Wit 
et al. 2009; Wheelwright et al. 2012). We decided to analyse 
images made from the nominal and orthogonal scans sepa¬ 
rately which maximises the resolution in the minor axis of 
the elongated PSF. 

We compared the Vesta PSF to that of the PSF ob¬ 
jects from the SDP field. This was to make sure that they 
were consistent with each other as the PSF changes because 
the angle between the scan and the detector axis (here¬ 
after array-to-map angle, a) changes with map direction 
(Lutz 2012). The values of the array-to-map angle used are 
a = +42.5° for the nominal direction, same as Vesta in fig- 
ure[l] and a = —42.5° for the orthogonal direction (Molinari 
et al. 2010a). Separate scan images of the PSF object V1362 
Aql are shown in figure [3] The structure in these images is 
similar to that in the Vesta image in figure \l\ Each Vesta 
image was rotated so that the major axis was horizontal and 
the scan direction pointing to the right, and then rebinned 
to the coarser pixel scale used in Hi-GAL as in figure [U A 
vertical slice was then taken along the minor axis, three pix¬ 
els wide and centred on the peak pixel. Each of the three 
columns were normalised to the central peak value and then 



RA (J2000) 


Figure 2. Image of the PSF star V1362 Aql in the naive map 
of the 1=30 region. Note the offset cross shape caused by the 
addition of the nominal and orthogonal scans each of which has 
an elongated PSF. 


a mean and standard deviation were taken from the three 
values at each offset position. 

The same procedure was applied to the nominal and 
orthogonal PSF objects after first subtracting a mean back¬ 
ground level determined in an annulus surrounding the ob¬ 
ject, and the intensity profile of the slices compared. The 
images were rotated so the positive scan leg£] point in the 
same direction as the rebinned Vesta image. To ensure the 
best comparison between objects at the subpixel level, a 2D 
Gaussian was fitted to the image and the slice was shifted 
so that the zero offset coincides with the Gaussian peak. 

Figure [5] shows that the intensity slices of the rebinned 
Vesta and the PSF stars V1362 Aql and IRAS 18491-0207 
drop to about 1% of the peak or about 1A" from the centre, 
which is where the uncertainties on the background level of 
the PSF objects become significant. Hence, since most of the 
MYSO targets are of similar or brighter flux than these PSF 
objects, comparing the MYSOs to these PSF objects would 
introduce more noise especially in the wings. Therefore, we 
compared our MYSOs to the much higher signal-to-noise 
image of Vesta from figure [T] 

It is also worth noticing that the minor axis of the PSF 
objects is independent of whether the object was present in 
one (e.g. IRAS 18491-0207 orthogonal) or both (e.g. IRAS 
18491-0207 nominal) scan legs. However, for the purpose of 
modelling the continuum emission (see Section 2} a Vesta 
image averaged with its reflection along the major axis was 
used for sources present in both scan legs. The difference 
along the minor axis slice is less than 10% between this 
averaged Vesta slices and the original one. 


1 For each map direction the telescope sweeps the sky in a pattern 
composed of several parallel scan legs, thus a leg can be either 
positive or negative depending on the direction with respect to 
the first scan leg. 
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Figure 3. Nominal (left) and orthogonal (right) scan image of the PSF star V1362 Aql. The flux scale is the same for both images. 


[h!] 



Figure 4. Image of Vesta (a = 42.5°) rotated and rebinned to 
Hi-GAL resolution showing the slice along the minor axis used in 
this work (red box). Each pixel has a size of 3.2" and the image 
is 29 x 23 pixels in size. 


pled objects, between 2 to 4 point sources were detected at 
8/im in the GLIMPSE/IRAC observations within the Her- 
schel 70/im resolution (~ 6 /x ), but > 50% of the emission is 
dominated by the MYSO at 8/xm and totally dominates in 
MIPS 24/im images. Therefore, in what follows we will not 
consider multiplicity as major concern. 

Comparing G030.8185+00.2729 in figure [6] to the naive 
map of the PSF star V1362 Aql in figure [2] clearly shows 
that the former is more extended compared to the PSF star. 
Similarly, figure [7| shows the nominal and orthogonal scan 
images for G030.8185+00.2729 and these clearly show ex¬ 
tended emission along the minor axis direction compared to 
the PSF star in figure [3j Intensity profiles for slices along 
the minor axis were constructed for each of the MYSOs as 
described above for the PSF objects. In figures l8l to flOl these 
slices (square blue points with error bars) are compared with 
those for Vesta (dashed cyan line). Again the 70/im emission 
from the MYSO is clearly more extended in the minor axis 
direction than the Herschel PSF for all the MYSOs in both 
scan directions. 


3 OBSERVATIONS 

3.1 70/im Imaging of RMS MYSOs 

There are a total of 12 RMS MYSOs in the Z=30° SDP 
field and 7 in the Z=59° field. We visually inspected each 
of these and found that only one in the Z=30° and two in 
the Z=59° field were sufficiently isolated from neighbour¬ 
ing sources and/or complex background (e.g. filamentary 
structures) emission to allow an investigation of their ex¬ 
tended emission. The parameters of the these MYSOs are 
given in table O Figure [6] shows the naive 70/im map of the 
brightest of the MYSOs, G030.8185+00.2729, and an ex¬ 
ample MYSO ignored for its complex background and the 
presence of a nearby object, G030.4117-00.2277. In our sam- 


3.2 Sub-mm Radial Profiles 

Although Herschel images of the same fields at 160, 250, 
350 and 500/xm are available, we limited the use of these 
images only for the SED. However, ground-based sub-mm 
observations are available which are at higher resolution. Az- 
imuthally averaged radial profiles at 870/rm were obtained 
from ATLASGAL images (Contreras et al. 2013) and at 
450/xm from the SCUBA Legacy Catalogues (Di Francesco 
et al. 2008) and are shown in figure fill Annulii of 1.5 times 
the pixel size (6" at 870/rm and 3" at 450/im) were chosen 
to compute the average flux in each bin of angular distance 
from the peak. The errors in each bin were estimated by the 
standard deviation of the fluxes in each annulus. 
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Table 2. Parameters of the isolated RMS MYSOs found within the two Hi-GAL SDP fields at 70/rni. 


Name 

RA 

Dec 

d 

(kpc) 

L 

(L©) 

Vh/d b 

(L© 1/2 kpc -1 ) 

./Vo/zin 

(Jy) 

/ 170//111 

(Jy) 

/ 250 /im 

(Jy) 

/ 350 /xm 

(Jy) 

/ 500 /^m 

(Jy) 

G030.8185+00.2729 

18:46:36.6 

-01:45:22 

5.7 

l.lxlO 4 

18.4 

321 

269 

131 

57.6 

22.7 

G058.7087+00.6607 

19:38:36.8 

+23:05:43 

4.4 

4.4xl0 3 

15.1 

30.9 

70.5 

44.4 

23.2 

12.2 

G059.8329+00.6729 

19:40:59.3 

+24:04:44 

4.2 

1.9xlO Sa 

10.4 

150 

361 

150 

61.0 

25.2 


a Note this object is in a cluster with several other YSOs within about 5". Its GLIMPSE 8p,m flux is only 20% of the larger MSX beam 
8pm flux and its total luminosity has therefore been reduced by this amount to reflect the fact there may be other luminosity sources in 
the large beam far-IR measurements of bolometric luminosity (see Mottram et al. 2011). 

b The physical size of a spherical dusty cloud heated to a particular temperature by a central source depends on the square root of the 
heating source luminosity which determines the spatial scales of the solution to the radiative transfer equation (Ivecic & Elitzur 1997). 
The angular size is then inversely proportional to the distance. Therefore, it is an indicator of how resolved is a source (see Section 0. 




Figure 5. A comparison of the intensity profile of Vesta (solid black line) and PSF stars V1362 Aql (dash-dotted red line) and IRAS 
18491—0207 (dashed green line) slices in the nominal (left) and orthogonal (right) map directions. Note the agreement within the errors 
between Vesta and the PSF objects out to about the 1% level. 
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Figure 6. Image of the RMS MYSOs G030.8185+00.2729 (left) and G030.4117-00.2277 (right) from the naive map of the 1=30 region. 
Note the partially resolved nature of G030.8185+00.2729 compared to the PSF star in figure [2] and the background level and morphology 
compared to G030.4117-00.2277. 
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Figure 7. Nominal (left) and orthogonal (right) scan images for the RMS MYSO G030.8185+00.2729. The flux scale is the same for 
both images. 


4 RESULTS 

To interpret the extended emission that we see at 70/im we 
have adapted the modelling procedure used by de Wit et al. 
(2009). We used the same grid of spherical radiative trans¬ 
fer models for MYSOs that was calculated by de Wit et al. 
(2009) using the DUSTY code (Ivezic & Elitzur 1997). The 
grid of 120 000 models spans a range in density law exponent 
p where n(r) oc r~ p with p varying from 0 to 2 in steps of 
0.5, Ay from 5 to 200 in steps of 5, and the ratio of outer ra¬ 
dius to sublimation radius, Y = Router/Rsub, varying from 
10 to 5000. For this study other model grid parameters were 
kept constant. These include the stellar effective tempera¬ 
ture that was kept at 25 000 K corresponding to a B0 V star 
as the IR emission is insensitive to this parameter. For the 
dust model we used the TSM’ model as described in de Wit 
et al. (2009) that consists of Draine & Lee (1984) graphite 
and silicate with an MRN size distribution (Mathis et al. 
1977). This has a sub-mm emissivity law with a slope of 
P = 2. The dust sublimation temperature was kept constant 
at 1500 K. 

Each model was scaled to the appropriate luminosity 
and distance for the MYSOs in table [2] A circular image 
of the emergent 70/im emission from the spherical model 
was generated and then convolved with the Vesta PSF re¬ 
binned for the Hi-GAL pixel scale. As before, an intensity 
profile slice was generated from an average of the three rows 
across the minor axis of the PSF direction normalised to 
the peak pixel. These model slices were then compared to 
the observed ones, both in the nominal and orthogonal scan 
directions. 

Simultaneously with fitting the intensity profile slice we 
also fitted the spectral energy distributions (SEDs). The lu¬ 
minosity used to scale the models comes from fitting the 
SED. The data points in the SEDs in figures [8] to [10] are 
from 2MASS (Skrutskie et al. 2006), GLIMPSE (Church- 
well et al. 2009), MSX (Price et al. 2001), Herschel (this 
work), sub-millimetre (Di Francesco et al. 2008; Contreras 


et al. 2013) and millimetre observations (Beuther et al. 2002; 
Beltran et al. 2006). Errors of 10% were adopted for all the 
SED data points to account for uncertainties in the abso¬ 
lute calibration across different datasets. During the fitting 
procedure reduced-x 2 (hereafter x 2 ) values were calculated 
for both the fits to the intensity slice and SED, where the 
degrees of freedom were the number of fitted points minus 
one. These are each placed in order of increasing x 2 and the 
model that is top of the combined order is considered to be 
the best fitting model (e.g. de Wit 2009). 

The best combined fits to the 70/xm intensity profile and 
SED for each direction are shown in figures [8] to [10] whilst 
the parameters are listed in table [3] In what follows, the re¬ 
sults are not referred to any particular scan direction unless 
otherwise stated. Reasonable combined fits to the intensity 
profile of most of the objects are obtained with the models, 
with x 2 near 1 for the 70/xm profile. The SED fit shows that 
fluxes at A < 3/im are always underestimated. This is com¬ 
mon for spherical models as they do not account for near-IR 
light being scattered and escaping from the bipolar outflow 
cavities (de Wit et al. 2010). The average power law index 
of the best fitting models is p — 0.5. 

Figure ITT1 show the profiles at 450 and 870/xm as seen 
in the combined fit of the SED and 70/xm intensity profile. 
Sub-mm radial profiles were also used instead of the 70/rni 
slices to analyse the effects of the spatial information from 
different wavelengths on the density distribution fit. These 
profiles constrain the density distribution of the cool outer 
regions. The combined 850/450/xm profile and SED best fits 
have an average exponent p = 0.5, which is consistent with 
the combined 70/im profile and SED fit. 

Finally, fits to each individual observation were also cal¬ 
culated. The results show that the best fits to radial profiles 
have exponents between 1 and 2, whilst the fits to SEDs have 
exponents between 0 and 0.5. In addition, the exponent of 
the best fitted models to the radial profiles is independent 
of Ay and the size of the cloud. 
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Figure 8. The combined best fit model (solid black) in terms of the 70p,m scan profile (top) and SED (bottom) compared to the data 
(blue squares) for G030.8185+00.2729. The left hand panels are for the scan in the nominal direction whilst the right hand panel is for 
the scan in the orthogonal direction. The Vesta PSF scan is shown in the top panel (dashed cyan) to illustrate the extended nature of 
the MYSO emission. An error of 10% of the total observed fluxes was considered in the observed SED. 


5 DISCUSSION 

The average power law index is shallower than the p ~ 1 
exponent in the sample of de Wit et al. (2009) who fitted 
24.5/im intensity profiles and SED. In addition, our aver¬ 
ages do not agree in general with other works which have 
found that the values of the power law index vary between 
1 and 2 by combining SED and sub-mm observations (e.g. 
Mueller et al. 2002). Nevertheless, in the particular case of 
G030.8185+00.2729, Williams et al. (2005) obtained an ex¬ 
ponent of 0.5 by using the SED and 850/xm radial profile, 
which agrees with our results. The results of Williams et al. 
(2005) where obtained by including in the SED points with 
A ^ 12/xm whilst Mueller et al. (2002) included points with 
A ^ 30/xm. We experimented with also only fitting data with 
A ^ 30/xmand found that in general values of exponents are 
0.5 lower than using the whole SED, and therefore the expo¬ 
nents are still between 0 and 1. We obviously do not expect 
to match the average results of Mueller et al. and de Wit et 
al. since our sample has only 3 objects. 

The power law indexes for the slice only and 870/xm only 
cases vary between 1.5 and 2, whilst in the 450/xm only cases 
its value is p — 1. These values are consistent with those 
obtained by Beuther et al. (2002), who found an average 
value of p — 1.6, even though they derived their values from 
a power law fitted to the radial profiles instead of doing the 


radiative transfer. Of course, these models do not fit the 
SED well. 

On the other hand, the power law indexes for the fits 
to the SED only vary between 0 and 0.5. This is similar to 
previous studies that use a dust emissivity law with a slope 
of /3 = 2 (e.g. Chini et al. 1986). As discussed by Hoare et 
al. (1991), a shallower dust emissivity law allows fits with a 
steeper density distribution. In fact, inspection of the SED 
fits at A > 100/xm in figure [9] appears to show that the A -2 
emissivity law used in the modelling is slightly too steep. 
Moreover, a study of the dust emissivity law in these two 
regions by Paradis et al. (2010) shows that the emissivity 
slope should be ~ 1.5 in l = 30° and ~ 1 in l — 59° for a dust 
temperature of 30 K. Our higher value of the slope would 
also explain the large values of Ay, for steeper emissivity 
laws need more dust mass to match the dust emission in the 
far-IR/submm. 

The values of Ay range between 95 and 200 for the 
combined SED and slice fit, but most of the sources have an 
Ay of 200. This is consistent with them being massive, young 
and embedded objects in their parental clouds. However, 
de Wit et al. (2009) found lower values than ours. In the 
particular case of G030.8185+00.2729, Williams et al. (2005) 
found a value ~ 4 times larger than ours, and using the 
method of Mueller et al. (2002) we obtain similar values of 
Ay as those obtained by considering all the points in the 
SED. Either way, all these results point towards values of 
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Figure 9. As in figure [8] but for G058.7087+00.6607. 


Ay greater than 90 magnitudes, and the inclusion of points 
at smaller wavelengths does not determine the value of Ay 
though it helps to constrain it. The value of Ay seems to be 
determined by the amount of dust necessary to reproduce 
the far-IR/sub-mm dust emission given its emissivity law. 

Table [2] shows the values of the y/L/d ratio, which has 
previously found to be a good indicator of how resolved these 
objects are (e.g. Wheelwright et al. 2011). This ratio is pro¬ 
portional to the angular size of the inner rim of the spherical 
envelope (Ivecic & Elitzur 1997) and, as is shown in table 02 
the envelopes have sizes a few thousands times the sublima¬ 
tion radius (~ 25 au for L = 10 4 L©), and should therefore 
be resolved at longer wavelengths. Figures 8 to 10 show the 
degree to which the objects are resolved at 70/im agrees with 
this. 


To explore whether this extends to the other 16 MYSOs 
with more complex background/neighbouring sources, we 
repeat the procedures used in the previous sections to ob¬ 
tain slices from the other sources in the Hi-GAL fields and 
from two models with similar physical properties as those 
obtained by the radiative transfer results, and then fitted ID 
Gaussian to measure the FWHM of these slices to see how 
resolved the sources are. Figure [12] shows the relation be¬ 
tween the FWHM and the y/h/d ratio. All observed sources 
are consistent with the models with some of them more ex¬ 
tended due to the complex background. 


6 CONCLUSIONS 

We have presented 70/im observations made with the Her- 
schel PACS instrument towards two regions of the Galactic 
plane and identified three relatively isolated MYSOs. The 
peculiarities of the Hi-GAL survey PSF and its effects on 
the MYSOs observations were analysed. The sources in our 
sample are all partially resolved at 70/xm. 

Using spherical radiative transfer models to simultane¬ 
ously fit the 70/rni profile and SED, we find we need a density 
law exponent of around 0.5. This is shallower than we previ¬ 
ously found from fitting partially resolved 24.5/im ground- 
based imaging, though both observations give an exponent 
between 0 and 1. It is also shallower than expected for in¬ 
falling material (p = 1.5). This could be due to rotational 
support, but since the emitting region is well outside of the 
expected disk/centrifugal radius (less than a few thousands 
AU, e.g. Zhang 2005) this is unlikely. It is more likely due 
to warm dust along the outflow cavity walls as seen in the 
mid-IR. We will investigate this further using 2D axisym- 
metric models. Intrinsic asymmetry could explain why we 
do not always get the same results on the same object from 
the nominal and orthogonal scan directions. 

Finally, the images at 70/im were smeared along the 
scan direction due to the scan speed. Moreover, the lack of 
PSF stars in the fields does not allow a characterisation of 
the PSF specific for these observations. Therefore, slow scan 
data, with a better behaved PSF, will be better to map and 
constrain the matter distribution of MYSOs. In particular, 
if the dust emission at 70/im comes from a non-spherical 
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Figure 10. As in figure [8] but for G059.8329+00.6729. 


Table 3. Best fit model parameters for the fits to the 70p,m intensity slice and the 450/rni and 870/xm radial profiles. 


Name 

Scan 

Fit 

V 

A v 

Y' a 

X.SED 


X450/an 

X870pm 

G030.8185+00.2729 

nominal 

SED + 70[im 

1.0 

200 

5000 

109 

1.4 

4.7 

0.23 



SED + 450/im 

0.5 

170 

2000 

40 

3.2 

0.4 

1.48 



SED + 870/un 

0.5 

120 

5000 

73 

3.6 

4.1 

0.04 


orthogonal 

SED + 70p,m 

0.5 

200 

1000 

77 

0.7 

1.9 

2.17 



SED + 450/un 

0.5 

170 

2000 

40 

1.6 

0.4 

1.48 



SED + 870/im 

0.5 

120 

5000 

73 

13.6 

4.1 

0.04 

G058.7087+00.6607 

nominal 

SED + 70fim 

1.0 

120 

5000 

132 

0.6 


0.48 



SED + 870/im 

0.5 

150 

5000 

62 

3.2 


0.25 


orthogonal 

SED + 70p,m 

0.0 

95 

5000 

45 

2.6 


0.35 



SED + 870/un 

0.5 

150 

5000 

62 

2.8 


0.25 

G059.8329+00.6729 

nominal 

SED + 70p,m 

0.5 

200 

1000 

4800 

4.8 


1.42 



SED + 870/un 

0.5 

200 

5000 

101 

6.2 


0.57 


orthogonal 

SED + 70p,m 

0.0 

200 

2000 

403 

1.5 


1.88 



SED + 870/im 

0.5 

200 

5000 

101 

2.0 


0.57 


Notes: The values of x 2 correspond to the reduced % 2 - 

a Y = Router / R S ub with Router the outer radius and R su b the sublimation radius of the envelope 


structure like bipolar cavity walls, data at 70/im can provide 
useful insights for 2D models. 
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Figure 12. Relation between y/L/d and FWHM of a ID Gaussian fitted to the 70/im slices from MYSOs in the l = 30° and l = 59° fields 
and models. The bar ranges are defined by the FWHM of the fit to the nominal and orthogonal directions, and the blue bars correspond 
to our sample. The predicted relation from two models with p = 0.5 and Ay = 200 is shown in red dashed line for Y = 1000 and in 
greed dot-dashed line for Y = 5000. Model images were convolved with the nominal Vesta PSF and a mean error of 1.2 /x is estimated 
for the Gaussian fit. The horizontal cyan dotted line represents the FWHM of the Vesta nominal slice. 
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